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Abstract. - We investigate a minimal model for non-crystalline water, defined on a Husimi 
lattice. The peculiar random-regular nature of the lattice is meant to account for the formation 
of a random 4-coordinated hydrogen-bond network. The model turns out to be consistent with 
most thermodynamic anomalies observed in liquid and supercooled-liquid water. Furthermore, 
the model exhibits two glassy phases with different densities, which can coexist at a first-order 
transition. The onset of a complex free-energy landscape, characterized by an exponentially large 
number of metastable minima, is pointed out by the cavity method, at the level of 1-step replica 
symmetry breaking. 



Introduction. — In the last decade, there have been 
' several attempts at describing structural glasses by means 
of simple lattice models, in the analytical framework of 
random-regular (Bethe or Husimi) lattices [iHl] • For these 
models, the cavity method generally predicts a discontin- 
uous replica-symmetry breaking, i.e., a sudden emergence 
of an exponentially large number of metastable free-energy 
minima. On the other hand, great interest has been at- 
tracted by polyamorphism {i.e., the existence of different 
glass forms of the same substance) [6j , both because of the 
relative novelty of the phenomenon (first discovered for 
water in 1985 [7^) and because of its relationship with the 
■ popular "second critical point" conjecture, put forward by 
\ Stanley and coworkers [8] to explain water anomalies [9]. 

Lattice models have long been used for investigating 
water [TUHH]. In some recent papers, coauthored by one 
of us [2T1[22], it has been shown that a first-order (quasi- 
chemical) approximation on a tetrahedral cluster is 
extremely effective to compute the phase diagram and the 
thermodynamic properties of a special class of water-like 
models, derived from the early Bell model [TU]. These 
models, defined on the (regular) body-centered cubic (bcc) 
lattice, do not contain a mechanism capable of inducing 
glassy behavior, so that they exhibit only crystalline (ice- 
like) phases at low temperature. 

The present letter is motivated by the following obser- 
vations, i) For a generic lattice model, the aforementioned 
quasi-chemical approximation, with a given choice of the 



associated cluster, coincides with the exact solution of a 
corresponding model, defined on a Husimi lattice made up 
of clusters of the same type p3ll24] , under the hypothesis of 
replica symmetry [251126] . As a consequence, we can define 
suitable Husimi lattice models, whose high-temperature 
behavior turns out to be very similar to that of the origi- 
nal water-like bcc-lattice models, ii) It is known from ex- 
periments that directional correlation of hydrogen bonds 
in real water is almost completely lost after the second 
consecutive bond (see [27] and references therein). The 
random nature of a Husimi lattice can roughly describe a 
similar scenario (this was not possible on a regular lattice), 
iii) We expect that the resulting frustration might hamper 
the onset of ice-like order, allowing for the possibility of 
glassy behavior at low temperature. Such a possibility can 
be easily investigated by the cavity method [25] . 

Concerning point i), let us remark the nontrivial dif- 
ference between a Husimi tree (which is, a system with a 
boundary) and a Husimi lattice (which is, a system with- 
out a boundary, in which all sites have the same coordi- 
nation number). The latter system locally exhibits the 
same tree-like structure as the former, but contains loops 
on a larger scale [25]. The two systems may be consid- 
ered equivalent as long as the frustration arising from the 
presence of loops is not strong enough to induce replica- 
symmetry breaking (i.e., glassy behavior). 

The main finding of this work is a Husimi lattice model 
predicting two different glassy phases, which we are led to 
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Fig. 1: Portion of the Husimi lattice. Each site is connected 
to 4 plaquettes. The bottom site, connected only to the cen- 
tral plaquette, gives a pictorial view of the "cavity bias" or 
"message" concept (see the text). 



identify with the low- and high-density amorphous (LDA, 
HDA) ices, observed experimentally [51 H]. This model 
provides interesting insights about the controversial na- 
ture of water polyamorphism and, more in general, the 
metastable phase diagram of water [9] . 

The model. — We consider a Husimi lattice made up 
of oriented plaquettes of 4 sites, as shown in fig. [T] Each 
site i (characterized by a configuration variable Xi) may 
be empty [xi — 0) or occupied by a water molecule, in 2 
possible configurations [xi = 1,2). Each water molecule 
possesses four equivalent bonding arms, which can point 
toward nearest neighbor sites, and may be either concor- 
dant {xi — 1) or discordant {xi = 2) with edge orienta- 
tions. An attractive potential energy — e < is assigned 
to any pair of nearest-neighbor occupied sites. This is the 
ordinary van der Waals contribution. A hydrogen bond 
is formed, yielding an extra energy —rj < 0, whenever 
two nearest-neighbor molecules point an arm toward each 
other, with no distinction between donors and acceptors 
(more precisely, a bond is formed on an edge oriented from 
i to j iS Xi — 1 and Xj — 2). The hydrogen bond may be 
weakened by the presence of extra molecules in the same 
plaquette, namely, each extra molecule gives rise to an en- 
ergy penalty rjc/2, with c G ]0, 1[. The hamiltonian of the 
system can be written as a sum over plaquettes 



(1) 



where site indices i,j,k,l are enumerated according to 
the plaquette orientation. The elementary contribution 
hx,y,z,w ("plaquette energy") can be written as 



"'X^y^z^w — ^x^y^z^w ~r ^y^z^w^x 

where 



"t~ '^z,w,x,y ~i~ '^w,x,y,z 5 (2) 
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Ttx ■ / -, f^z + n-w 1 /„N 

'fJ-— erixHy ~ Tjhx^y 1 - c ) , (3) 



nx is an "occupation function" {ux = Oifa:; = 0;nj, = l 
otherwise), bx.y is a "bond function" {bx.y = 1 if x = 1 and 



y ~ 2] bx.y — otherwise), and /i is the chemical poten- 
tial (we study the grand-canonical ensemble) . Looking at 
eq. ([2]), one immediately realizes that the plaquette energy 
is invariant under circular permutations of the configura- 
tion variables. As a consequence, the full hamiltonian ([T]) 
is imaffected by the choice of the "first site" i in each 
plaquette. 

For the sake of brevity, we have given only a formal 
description of the model. Indeed, the latter is the Husimi- 
lattice version of a water-like model similar to that pro- 
posed by Roberts and Debenedetti [13]. This can be eas- 
ily deduced by comparing the original paper [T3] with 
the tetrahedral cluster approximation developed in |18j . 
The square plaquettes of the Husimi lattice correspond to 
tetrahedral clusters on the bcc lattice, whereas plaquette 
orientations are related with the geometric structure of 
model molecules. Let us remark that, even though the 
current model neglects some details considered in |14| and 
[18] (namely, the distinction between donors and accep- 
tors [13], and the presence of extra nonbonding config- 
urations [131 [H]), it nonetheless reproduces the typical 
thermodynamic anomalies observed in real water. Such 
a result confirms that the physical mechanism underlying 
the anomalies is mainly based on the "weakening param- 
eter" c, which favors states characterized by a positive 
correlation between higher local entropy (weaker bonds) 
and higher local density [T311I3. A standard statistical- 
mechanical argument [5] relates such a positive correlation 
to a negative thermal expansion coefficient, i.e., to the on- 
set of a density maximum. 

The replica-symmetric (RS) solution. Accord- 
ing to the cavity method [53], in the RS assumption, the 
model can be solved by a recursion equation for so-called 
"cavity biases" or "messages". An elementary message 
to""*' represents the probability of the x configuration at 
the i site, when the latter is detached from all plaquettes 
except a. With reference to fig. [H denoting by frix the 
message from the central plaquette to the bottom site, 
the recursion equation reads 



EEE' 

a^i— 3^2—0 3^3—0 



i—l a—1 



(4) 



where /3 — l/k^T is the inverse temperature and / is a 
normalization constant, ensuring that ^^^qTtIx = 1- Let 
us note that, in our case, since the lattice has no local 
heterogeneities, the messages do not depend on the po- 
sition (one can drop the a ^ i superscript). Therefore, 
eq. Q eventually simplifies to a set of three equations 
{x — 0, 1, 2) of the fixed-point form 
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where we have omitted the trivial normalization constant. 
These equations can be solved numerically by simple iter- 
ation. The actual probability px of the x configuration can 



p-2 



A discrete model of water 



then be evaluated by considering the operation of attach- 
ing four equivalent branches, like that depicted in fig.[Il to 
the same "root" site [24]. One obtains Px oc nix'^, where 
we have again omitted a normalization constant, needed 
to ensure "^^^qPx = 1- The RS solution is always char- 
acterized by pi = P2, *-e-, no preference between the two 
molecule configurations. The mass density p is simply re- 
lated to the occupation probability Pocc = 1 — Po — P1+P2 
as p = (M/v) pocc, where M w 18g/mol is the molecular 
mass of water and v is the volume per site (which is an 
adjustable parameter of the model). 

The grand-canonical free energy per site F/N (let F 
and N denote respectively the thermodynamic free-energy 
times /3, and the number of lattice sites) can be evaluated 
as a function of the messages as [53] 
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(6) 

Pressure can then be computed as P = —{kQT/v){F/N). 
In the presence of multiple solutions of eq. ([S]), i.e., of 
competing phases, the thermodynamically stable one is se- 
lected by the lowest free-energy (highest pressure) value. 
In general, the knowledge of the cavity biases allows one 
to compute all the thermodynamic properties of interest, 
including response functions, such as the specific heat, the 
isothermal compressibility kt = {din p/dP)T, andtheiso- 
baric thermal expansion coefficient ap — —{din p/dT)p, 
whose anomalous behavior is specially relevant in real 
water. It is also possible to obtain numerically simple 
equations for the loci of divergence of the response func- 
tions (spinodals) and the temperature of maximum den- 
sity (TMD) locus at constant pressure, defined by ap = 0. 

A fitting procedure has been performed to fix the model 
parameters. The hydrogen bond energy 77 and the volume 
per site v are easily computed as a function of the van der 
Waals energy e and the weakening parameter c, by requir- 
ing the maximum density at 1 atm to be 1 g/cm^ and to oc- 
cur exactly at 3.984 °C [H]- We have then fitted only the 
two parameters e and c by minimizing (in a least-square 
sense j3 the difference between a set of experimental den- 
sity values at 1 atm (fig. [S] inset) and the corresponding 
theoretical values. The final result is: 77 w 7.786 kJ/mol, 
e « 3.621 kJ/mol, c « 0.6450, v w 15.70 cm^/mol. 

In fig. [5] we report the RS phase diagram in the 
temperature-density plane. We obtain a quite good agree- 
ment with experimental data for the whole liquid-vapor 
binodal curve (not fitted). The critical exponent is not 
correct, due to the mean-field nature of the model, but the 
critical density pc ~ 0.3417 g/cm'^ is remarkably close to 
the experimental value 0.322 g/cm^ [55]. The TMD locus 
correctly displays a negative slope, quantitatively similar 




^We have used the MATLAB routine Isqnonlin. 
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Fig. 2: Temperature-density phase diagram in the RS assump- 
tion. Thick solid lines denote binodals. Thin lines are defined 
by the following conventions: A solid line denotes the liquid 
spinodal, a dash-dotted line denotes the TMD locus, a dashed 
line denotes the stability limit of the RS solution, and a dotted 
line denotes the zero-entropy locus. Solid and open squares 
denote experimental results for the liquid-vapor binodal [29j 
and the TMD locus '30 (the latter data refer indeed to heavy 
water; pressure values have been converted to density values 
by [29]). The inset displays the isobar at P = 1 atm, along 
with the experimental values [28] used for parameter fitting. 



to the experimental one. Moreover, in the metastable re- 
gion of the liquid phase at very low densities (negative 
pressures), the TMD locus exhibits a slight reentrance 
{i.e., positive slope), which has been predicted by simu- 
lations with various intermolecular potentials [3TJ|ni]. At 
odd with such simulations, our TMD line eventually meets 
the spinodal line. In the pressure-temperature diagram, 
the two curves meet tangentially at a pressure minimum, 
as required for thermodynamic consistency [33 , but the 
spinodal is not found to re-enter the positive pressure re- 
gion. The latter fact contradicts Speedy's early conjec- 
ture [34] about the divergent-like behavior of response 
functions in the supercooled liquid regime. Conversely, 
our model seems to support Stanley's conjecture, as it pre- 
dicts, at very low temperature, a "second critical point" , 
terminating a coexistence region of two different liquid- 
like phases. Unfortunately, all this coexistence region lies 
below the stability limit of the RS solution (dashed line 
in fig. [2]), where the latter is no longer valid. The inad- 
equacy of such a solution in this regime is also pointed 
out by the fact that its entropy becomes negative below a 
given temperature (dotted line in fig. |2]). 

Fig. 13] displays the isobaric specific heat and the isother- 
mal compressibility as a function of temperature, at at- 
mospheric pressure. The well-known anomalous behav- 
ior of these response functions, related respectively to en- 
tropy and density fluctuations [Qi, is reproduced by the 
model in a qualitatively correct fashion. This result is 
achieved without the singularity required by Speedy's con- 
jecture. Indeed, upon decreasing temperature, the the- 
oretical curves do not exhibit any divergence, but only 
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Fig. 3: Normalized isobaric specific heat (left) and isother- 
mal compressibility (right) as a function of temperature at 
P = 1 atm. Lines and symbols denote respectively theoretical 
and experimental 28; results. The theoretical (experimental) 
minima occur at 37.45 °C (36 °C) and 85.58 °C (46.5 °C). 



pronounced maxima. Quantitative agreement with exper- 
iments is not so good in this case, so that we report data 
normalized to the minimum values c^'"'' and k^™'*. 

The glassy phases. — Let us now investigate what 
happens when the RS ansatz does not hold. We make 
use of the cavity method at the level equivalent to 1-step 
replica symmetry breaking (RSB) |25| . For a given free- 
energy landscape, the complexity function S(i^) is defined 
as the log- number of minima ( "pure states" , or simply 
"states" ) with free energy in the range [F, F + dF] . If 
the number of states is exponentially large, then S is an 
extensive quantity. One can define a pseudo free-energy 
("repHcated" free-energy) as [35] 



(7) 



where Fa is the free energy of the a state, the sum runs 
over all states, and ip is a pseudo inverse-temperature 
(Parisi parameter). A saddle-point calculation shows that 
$(■0) is the Legendre transform of S(F), so that $(V') 
allows one to reconstruct S(-F) via the parametric repre- 
sentation 
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The replicated free-energy (and all the thermodynamic 
quantities) can be computed, as a function of '0, by solving 
an integral equation for the message distributiord V{m) 
over the states 
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(9) 

where both m and / are defined by eq. ^ as functions of 
the set of incoming messages {m°'^^}, and (f> ensures the 
normalization condition J 'P{m)dm — 1. Equation ^ 
can be solved numerically by a population dynamics tech- 
nique |25j . 



^Hereafter, m without subscript denotes a normalized 3- 
component array (mo,mi,m2). 




Fig. 4: F/N (dashed line) and {F - E) /N (solid line) as a 
function oiip &tT = 135 K, yi/rj — —2.38. Thinner lines denote 
unphysical parts of the curves. 

In the RS phase, V{m) is a delta-function and eq. (|4]) is 
recovered. The stability limit of the RS phase, displayed 
in fig. [21 can be determined by assuming that V{m) has 
only a slight variance around the RS value. One obtains a 
linearized "propagation equation" for the covariance ma- 
trix of V{m) 



x'=0 y'=0 



(10) 



where the (9 x 9) "transfer matrix" is 

drhx 



K. 
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x,y;x\y' 
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(11) 



and the Jacobians drhx/dm'^7^^ can be evaluated from 
eq. ([4|) . The RS solution becomes unstable when the maxi- 
mum eigenvalue of the transfer matrix (|11|) becomes larger 
than 1. 

In the RSB phase, the thermodynamically relevant 
states are identified by a minimum (wrt "0) of the func- 
tion F — Yj^ which properly takes into account that the 
probability measure may be split over a large number of 
states with free energy around F. We report the typi- 
cal behavior of F/N and {F - E) /N in fig. [H The left- 
most part of the curves {ip < "0*) is unphysical, because 
dF/d'0 — d^^/dV'^ < 0, which is inconsistent with eq. ([7|). 
The rightmost part (0 > ipo) is also unphysical, because 
S < 0. Therefore, the physical minimum of — E is at- 
tained at the crossing point 0o, such that E('0o) =0. In 
the jargon of replica theory, this is called a "condensed" 
glass phase. The number of thermodynamically relevant 
states is sub-exponential, but there exists an exponen- 
tially large number of metastable states. Fig. [5] shows 
that the former states exhibit the largest density, p{ipo)y 
whereas the latter span a range of smaller density values, 

P(0*) < P< P(0o)- 

Upon increasing temperature, we observe that the RSB 
phase undergoes a continuous transition to the RS phase. 
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Fig. 5: Parametric plot of T,{ip)/N vs. p{ip) at T = 135 K, 
fi/r) = —2.38 (as in fig. 2)). Thinner lines denote unphysical 
parts of the curve. 

at odd with the glass models cited above [l]-[5]. As pre- 
viously mentioned, those models predict (upon decreasing 
temperature) an abrupt appearance of many metastable 
states ("dynamical transition"), where the system is ex- 
pected to get stuck during quenching. Such a difference 
is rather intriguing. In real water, the glass transition ap- 
pears to be much "weaker" than that found in the usual 
molecular liquids [55], and the very existence of a glass 
transition in water (before recrystallization upon heating) 
has been much debated _l37j[38'. This "weakness" is be- 
lieved to reflect a greater ease for the system to explore 
its energy landscape, i.e., a greater accessibility of the 
ideal glass state Qualitatively speaking, such a sce- 
nario might be consistent with the absence of a dynamical 
transition, as predicted by our model. 

Let us also consider the RSB phase at constant temper- 
ature. Fig.|B]shows that, upon increasing the chemical po- 
tential, a higher density solution appears discontinuously. 
It turns out that even the latter is a condensed glass phase, 
in the previously explained sense. As previously men- 
tioned, we find it natural to identify these phases with the 
LDA and HDA ices. Remarkably, the density range pre- 
dicted by the RSB solution for LDA is consistent with ex- 
perimental values, at odd with the low-density (unstable) 
RS solution. The behavior of HDA, as a function of the 
Parisi parameter if), is qualitatively similar to that of LDA, 
displayed in fig. ID Nevertheless, the density values of the 
metastable states are higher than the density of the ther- 
modynamically relevant states, i.e., p{ip*) > p{i'o) (the 
p-E plot is a mirror image of that of fig. [5]) . Such a result 
guarantees that, even in the region where both solutions 
exist, they are always separated by a density gap with no 
states. As a consequence, our model supports the picture 
of a real first-order-like LDA-HDA transition [6, 9j, ex- 
cluding the possibility of a smooth crossover, suggested by 
some experimental observations |39| . The transition pres- 
sure at T = 135 K turns out to be about 28.6 MPa, much 
lower than the experimental estimate (^ 200 MPa) [lUj . 
but the decreasing trend, upon increasing temperature, is 



Fig. 6: Density p as a function of the normahzed chemical 
potential /i/r; at T = 135 K. Thinner hues denote the (unsta- 
ble) RS solution; thicker lines denote p{ipo) {i.e., the density 
of the thermodinamically relevant states) in the RSB solution. 
Dotted lines denote metastable continuations. Experimental 
density ranges of LDA and HDA ices [40] are highlighted on 
the left axis. 

correctly reproduced. 

The transition line terminates in a critical point, which 
occurs around T « 198 K, P ^ 16.7 MPa, and above which 
the two glassy phases become indistinguishable. Note 
that, since LDA and HDA are, respectively, characterized 
by p{ip*) < p(V'o) and p{ip*) > p{ipo), indistinguishability 
implies p{ip*) = p(V'o) in the supercritical region. This 
fact gives rise to a third kind of RSB phase (A), char- 
acterized hy ipo — ip* = (see footnotqfl). The maximal 
complexity S(V'*) of LDA and HDA turns out to vanish in 
a continuous way, defining two more (LDA-A and HDA- A) 
transition lines. Fig.[7]shows numerical evidence that such 
lines meet precisely at the critical point. The A-phase may 
still be considered a glassy phase, since the message distri- 
bution 'P(m) is still nontrivial. As previously mentioned, 
the vanishing complexity means that such a distribution 
represents a sub-exponential number of pure states. As 
shown in fig. [71 the A-phase region extends until the RS 
stability limit, where V(m) eventually degenerates into a 
Dirac delta, and the RS solution becomes stable. 

Summary and conclusions. We have studied a 
simplified model of water, defined on a Husimi lattice, 
which can be solved by the cavity method. Despite its 
simplicity, the model seems to capture a lot of features 
arising in real (and simulated) water, including thermody- 
namic anomalies of the liquid phase and a low-temperature 
glassy phase displaying polyamorphism. To the best of 
our knowledge, this is the first lattice model describing all 
these features together. 

The glassy behavior originates in the competition be- 
tween the tendency of the model to form regular hydrogen- 
bond networks, driven by the energy term —rjb^ y of 
eq. ((31), and the topological disorder of the Husimi lat- 

^The corresponding complexity is T, = 0; ip values other than 
result in a negative complexity. 
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Fig. 7: RSB phase diagram (see the text). The dashed line 
denotes the RS stabihty hmit. Other hnes are 2nd order poly- 
nomial interpolations of symbols. 



tice, which induces frustration. In other words, the lattice 
randomness plays a relevant role in the model behavior, 
as it is responsible for suppressing the ice-like phases that 
appear in the corresponding regular lattice models. Even 
though the latter feature might be considered a weak- 
ness of the model, we believe that it roughly illustrates 
the physical mechanism underlying glassy behavior in a 
hydrogen-bonded (network-forming) fluid. In this frame- 
work, the onset of polyamorphism is related to the fact 
that, due to the weakening term +'i]bx^yc{nz + nju)/2, the 
original model predicts a low density network structure 
(stable at low pressure), besides the "ordinary" crystal- 
like ground state. The glass-glass transition may thus be 
viewed as a "frustrated" reminiscence of the coexistence 
between these two phases. 

In the future, it might be interesting to investigate 
whether a similar mechanism could give rise to glass-glass 
transitions even for lattice models predicting a dynamical 
transition scenario [THl]. A more technical issue, still de- 
serving investigation, is the stability of the current picture 
with respect to further RSB steps. 
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